% By Martin Andreasen, 2009
% This function computes bond prices up to a given maturity, provided the
% first bond prices is reported in the function g(x,sig) and stored at the
% position ny_1bond. This is the special version were a log-transformation
% is used. This implementation allows for no-zero third moment of
% the shock distribution.
% In this version, symmetries in Pxx and Pxxx are used

function [P,Px,Pxx,Pxxx,Pss,Pssx,Psss] = Get_Bond_Prices_Mom3_Log(M,Mxp,Myp,...
    Mxpx,Mxpxp,Mxpy,Mxpyp,Mypx,Mypxp,Mypy,Mypyp,...
    gx,gxx,gxxx,gss,gssx,gsss,hx,hxx,hxxx,hss,hssx,hsss,eta,vectorMom3,Maturity,ny_1bond,order_app,solveRiskTerms)


% The arrays
[ny,nx] = size(gx);
ne      = size(eta,2);
P       = zeros(Maturity,1);
Px      = zeros(Maturity,nx);
Pxx     = zeros(Maturity,nx,nx);
Pxxx    = zeros(Maturity,nx,nx,nx);
Pss     = zeros(Maturity,1);
Pssx    = zeros(Maturity,nx);
Psss    = zeros(Maturity,1);

% The all the derivatives for the first bond price
P(1,1)       = M;
Px(1,:)      = gx(ny_1bond,:);
Pxx(1,:,:)   = squeeze(gxx(ny_1bond,:,:));
Pxxx(1,:,:,:)= squeeze(gxxx(ny_1bond,:,:,:));
Pss(1,1)     = gss(ny_1bond,1);
Pssx(1,:)    = gssx(ny_1bond,:);
Psss(1,1)    = gsss(ny_1bond,1);

% The value of bond prices in deterministic steady state steady state
k = 2:Maturity;
P(k,1) = M.^k;


% The first order terms
% Recall that R(P(0,1)) = 1
for k=2:Maturity
    % With matrix notation
    Px(k,:) = Px(1,:)+ Px(k-1,:)*hx;
end
% Stop if order_app == 1
if order_app == 1
    return;
end


% The second order terms, Pxx
for k=2:Maturity
    tmp = zeros(nx,nx);
    for gama1 = 1:nx
        tmp = tmp + Px(k-1,gama1)*squeeze(hxx(gama1,:,:));
    end
    Pxx(k,:,:) = squeeze(Pxx(1,:,:)) + hx'*squeeze(Pxx(k-1,:,:))*hx + tmp;
end
% The constant correction, Pss
for k=2:Maturity
    Pss(k,1) = Pss(1,1) + Pss(k-1,1) + Px(k-1,:)*hss + trace(eta'*squeeze(Pxx(k-1,:,:))*eta) ...
        + Px(k-1,:)*eta*eta'*Px(k-1,:)' + 2*M^-1*Myp(1,:)*gx*eta*eta'*Px(k-1,:)'...
        + 2*M^-1*Mxp(1,:)*eta*eta'*Px(k-1,:)';
    
end
% Stop if order_app == 2
if order_app == 2
    return;
end


% The third order order terms, Pxxx
for k=2:Maturity
    for alfa1=1:nx
        for alfa2=alfa1:nx
            for alfa3=alfa2:nx
                tmp = 0;
                for gama1=1:nx
                    tmp = tmp + hx(gama1,alfa1)*hx(:,alfa2)'*(squeeze(Pxxx(k-1,gama1,:,:)))*hx(:,alfa3);
                end
                Pxxx(k,alfa1,alfa2,alfa3) = Pxxx(1,alfa1,alfa2,alfa3) + tmp ...
                    + hx(:,alfa1)'*squeeze(Pxx(k-1,:,:))*squeeze(hxx(:,alfa2,alfa3)) ...
                    + reshape(hxx(:,alfa1,alfa3),1,nx)*squeeze(Pxx(k-1,:,:))*hx(:,alfa2) ...
                    + reshape(hxx(:,alfa1,alfa2),1,nx)*squeeze(Pxx(k-1,:,:))*hx(:,alfa3) ...
                    + Px(k-1,:)*squeeze(hxxx(:,alfa1,alfa2,alfa3));
                
                % Using symmetry for alfa1 and alfa2
                if alfa1 == alfa2 && alfa2 ~= alfa3 %alfa1==alfa2~=alfa3
                    %Pxxx(k,alfa1,alfa1,alfa3)= Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa1,alfa3,alfa1) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa3,alfa1,alfa1) = Pxxx(k,alfa1,alfa2,alfa3);
                end
                % Using symmetry for alfa2 and alfa3
                if alfa1 ~= alfa2 && alfa2 == alfa3  %alfa1~=alfa2==alfa3
                    %Pxxx(k,alfa1,alfa2,alfa2)= Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa2,alfa1,alfa2) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa2,alfa2,alfa1) = Pxxx(k,alfa1,alfa2,alfa3);
                end
                % Using symmetry for alfa1,alfa2, and alfa3
                if alfa1 ~= alfa2 && alfa1 ~= alfa3 &&  alfa2 ~= alfa3 %alfa1~=alfa2~=alfa3
                    %Pxxx(k,alfa1,alfa2,alfa3) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa1,alfa3,alfa2) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa3,alfa1,alfa2) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa3,alfa2,alfa1) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa2,alfa3,alfa1) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa2,alfa1,alfa3) = Pxxx(k,alfa1,alfa2,alfa3);
                end
            end
        end
    end
end

if solveRiskTerms == 1
    % The third order order terms, Pssx
    % First some auxiliary variables
    Myp_gxx = zeros(nx,nx);
    for gama1=1:nx
        for gama2=1:nx
            Myp_gxx(gama1,gama2) = Myp(1,:)*squeeze(gxx(:,gama1,gama2));
        end
    end
    for k=2:Maturity
        tmp1 = zeros(1,nx);
        for beta1=1:ny
            tmp1 = tmp1 + 2*M^-1*Myp(1,beta1)*Px(k-1,:)*eta*eta'*squeeze(gxx(beta1,:,:))*hx;
        end
        tmp2 = zeros(1,nx);
        for gama1=1:nx
            tmp2 = tmp2 + eta(gama1,:)*eta'*squeeze(Pxxx(k-1,gama1,:,:))*hx;
        end
        Pssx(k,:) = -2*M^-1*Myp*gx*eta*eta'*Px(k-1,:)'*Px(1,:) ...                                       %1)
            -2*M^-1*Mxp*eta*eta'*Px(k-1,:)'*Px(1,:) ...                                          %2)
            +Pssx(1,:) ...                                                                       %3)
            +2*M^-1*Px(k-1,:)*eta*eta'*gx'...                                                    %4)
            *(Mypyp*gx*hx + Mypy*gx + Mypxp*hx + Mypx) ...
            +tmp1 ...                                                                            %5)
            +2*M^-1*Px(k-1,:)*eta*eta'*(Mxpyp*gx*hx+Mxpy*gx+Mxpxp*hx+Mxpx) ...                   %6)
            +2*M^-1*Myp*gx*eta*eta'*squeeze(Pxx(k-1,:,:))*hx ...                                 %7)
            +2*M^-1*Mxp*eta*eta'*squeeze(Pxx(k-1,:,:))*hx ...                                    %8)
            +Px(k-1,:)*eta*eta'*squeeze(Pxx(k-1,:,:))*hx ...                                     %9)
            +Px(k-1,:)*eta*eta'*squeeze(Pxx(k-1,:,:))*hx ...                                     %10)
            +tmp2 ...                                                                            %11)
            +hss'*squeeze(Pxx(k-1,:,:))*hx ...                                                   %12)
            +Px(k-1,:)*hssx ...                                                                  %13)
            +Pssx(k-1,:)*hx;                                                                     %14)
    end
    
    % The third order order terms, Psss
    for k=2:Maturity
        % Term b1)
        Psss(k,1) = Psss(1,1);                                                                                        %b1)
        
        % The terms b2), b3), b5), b6), b7), b8) and b9)
        tmp = 0;
        for phi1=1:ne
            tmp = tmp + 3*M^-1*(gx*eta(:,phi1)*vectorMom3(phi1))'*Mypyp*gx*eta(:,phi1)*Px(k-1,:)*eta(:,phi1)...            %b2)
                + 6*M^-1*(gx*eta(:,phi1)*vectorMom3(phi1))'*Mypxp*eta(:,phi1)*Px(k-1,:)*eta(:,phi1)...               %b3)
                + 3*M^-1*(eta(:,phi1)*vectorMom3(phi1))'*Mxpxp*eta(:,phi1)*Px(k-1,:)*eta(:,phi1)...                  %b5)
                + 3*M^-1*(Myp(1,:)*gx*eta(:,phi1) + Mxp(1,:)*eta(:,phi1))*vectorMom3(phi1)*Px(k-1,:)*eta(:,phi1)...  %b6)
                *Px(k-1,:)*eta(:,phi1)...
                + 3*M^-1*(Myp(1,:)*gx*eta(:,phi1) + Mxp(1,:)*eta(:,phi1))*vectorMom3(phi1)...                        %b7)
                *eta(:,phi1)'*squeeze(Pxx(k-1,:,:))*eta(:,phi1)...
                + Px(k-1,:)*eta(:,phi1)*Px(k-1,:)*eta(:,phi1)...                                                     %b8)
                *Px(k-1,:)*eta(:,phi1)*vectorMom3(phi1)...
                + 3*Px(k-1,:)*eta(:,phi1)...                                                                         %b9)
                *(eta(:,phi1)*vectorMom3(phi1))'*squeeze(Pxx(k-1,:,:))*eta(:,phi1);
        end
        Psss(k,1) = Psss(k,1) + tmp;
        
        % Term b4)
        tmp = 0;
        for beta1=1:ny
            for phi1=1:ne
                tmp = tmp + 3*M^-1*(eta(:,phi1)*vectorMom3(phi1))'*Myp(1,beta1)*squeeze(gxx(beta1,:,:))*eta(:,phi1)...
                    *Px(k-1,:)*eta(:,phi1);
            end
        end
        Psss(k,1) = Psss(k,1) + tmp;
        
        % The term 10b)
        tmp = 0;
        for gama1=1:nx
            for phi1=1:ne
                tmp = tmp + eta(:,phi1)'*squeeze(Pxxx(k-1,gama1,:,:))*eta(:,phi1)*eta(gama1,phi1)*vectorMom3(phi1);
            end
        end
        Psss(k,1) = Psss(k,1) + tmp;
        
        % Terms b11 and b12
        tmp =  Px(k-1,:)*hsss(:,1) ...                                                              %b11)
            + Psss(k-1,1);                                                                        %b12)
        Psss(k,1) = Psss(k,1) + tmp;
    end
end